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1 Introduction 



> : 

Lattice QCD was invented thirty years ago but only in the last few years has it finally ful- 
■ filled its promise as a precision tool for calculations in hadron physics. This review will 
Q> ' cover the fundamentals of discretising QCD onto a space-time lattice and how to reduce 
the errors associated with the discretisation. This 'improvement' is the key that has made 
the enormous computational task of a lattice QCD calculation tractable and enabled us 
to reach the recent milestone of precision calculations of simple 'gold-plated' hadron 
| masses. Accurate decay matrix elements, such as those for leptonic and semileptonic de- 
'"" j" 1 • cay of heavy mesons needed by the B factory experimental programme, are now within 
sight. I will describe what goes into such calculations and what the future prospects and 
D • limitations are. 



in 
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^ ! 2 Lattice QCD formalism and methods 
2.1 The Path Integral 

Lattice QCD is based on the path integral formalism. We can demonstrate this formalism 
by discussing the solution of the quantum mechanical problem of a particle moving in 
one dimension (Lepage 1998a). This could be solved using Schrodinger's equation with 
H = p 2 /2m + V(x) and [x,p] — i (in 'particle physics units', h = c = 1). We can also 
solve it using the path integral formulation and this is the basis for lattice QCD. 

A key quantity is the transition amplitude between eigenstates of position at, say, 
time t = U and tf. This is expressed as a functional integral over all possible paths x(t) 
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from t = ti to tf weighted by the exponential of the action, S, which is the integral of 
the Lagrangian, L. 



(x f {t } )\ Xl (t t )) = / vx{ty s ^ 



S[n 



dtL(x, x) — I dt[ 



,mx(t) 2 



V(x(t))}. 



(1) 

(2) 



The path integral can be evaluated by discretising time into a set of points, tj = ti + ja 
for j — 0,1 ... N and a the lattice spacing = (tf — ti)/N. The path x(t) then becomes 
a set of variables, Xj, and Equation^above requires us to integrate over each one, ie the 
problem reduces to an ordinary integral but over an N — 1-dimensional space. The end 
points, xo and a; at are kept fixed in this example. The most likely path is the classical one, 
which minimises the action (mx = V'), but the path integral allows quantum fluctuations 
about this, see Figure^ 




Figure 1. Discretised classical and possible quantum mechanical paths from Xi to Xf 
for a particle moving in one dimension. 

The i in front of the action in the exponential gives rise to the problem, both concep- 
tual and numerical, of adding oscillating quantities together. It is simpler to rotate the 
time axis to Euclidean time, t — > —it. Then the path integral becomes 



A 



Vx(t)e~ s[x 

S[x] = / L(x,x) = ^2 

Jti __n 



dx\dx2 ■ ■ ■ c?xjv-ie 
N-l 



-S[x] 



mX : 



V(xj)], 



(3) 



3=0 



where the integrals over the intermediate points in the path, x(t), are now explicit. We 
will ignore the normalisation of the integral, A. To perform the integral we must discre- 
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tise the Lagrangian so that it takes values at the discrete time points. Then 

N-l 

S=£[^(*i + i-Zi) 2 + an^)]- (4) 

3=0 a 

Clearly the accuracy of our discretisation will depend on a being small. However, as a 
becomes smaller at fixed physical time length, the number of points, N, increases and so 
does the computational cost. 

For large N an efficient way to perform the integral is by Monte Carlo. A set of 
possible values for Xj, j — 1, N is called a configuration. The configurations with most 
weight in the integral are those with large e . For maximum efficiency we want to 
generate configurations with probability e~ s - this is known as importance sampling. 
A simple method for doing this is the Metropolis algorithm. This starts with an ini- 
tial configuration {eg Xj all zero or chosen randomly). It then passes through the Xj in 
turn proposing a change to a given Xj of size e, ie a random number between — e and 
e temporarily added to Xk, say. The change in the action as a result of Xk changing is 
calculated. Call this AS. Note that this calculation only involves the Xj in the neighbour- 
hood of Xk and connected to it through the action. If AS < the change is accepted. 
If AS* > another random number, uniformly distributed between and 1 is generated. 
If e" AS > rand then the change is accepted. If not, Xk reverts to its previous value. 
(Lepage 1998a, di Pierro 2000) 

In this way a new configuration is generated after each sweep through the lattice. A 
set of configurations is called an ensemble. Calculations of various functions of the Xj on 
an ensemble then yield the physics results that we are after, such as the quantised energy 
levels available to the particle. 

In wavefunction language we can express 

(x f \e- iH ^-^\ Xi ) =J2r n (xf)M^)e- iEnitf - U) , (5) 

a 

by inserting a complete set of eigenstates of the Hamiltonian, H. When rotated to Eu- 
clidean time, and taking tf — U — T and x, = Xf = x we have 

(x\e- HT \x) = r„(x)Mx)e- EnT , (6) 

n 

where ip n (x) — (n\x). It now becomes clear that the result will be dominated by the 
ground state as T becomes large, because all higher states are exponentially suppressed. 
Integrating over the initial and final x values gives: 

dx(x\e- HT \x) -> e - EoT ,T^ oo. (7) 

This would be the result, up to a normalisation, of integrating by Monte Carlo over all 
the xj, setting the initial and final values to be the same, ie using periodic boundary 
conditions. The result looks very similar to a problem in statistical mechanics when 
formulated in this way, as did Equation [3] and this is not an accident. Indeed we can 
work at non-zero temperature if we take T finite, but here we will concentrate on the 
zero temperature case where in principle T — > oo, and in practice is large. 
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To investigate excitations above the ground state, we must interrupt the propagation 
of the ground state by introducing new operators at intermediate times. For example, 

(x(T)\x{t 2 )x(ti)\x(0)) = J Vx x(t 2 )x(h)e- s W 
(x(T)\x(0)) ~ fVxe- s W 

= KBbN^^e-^-^^-^.ta-ti-cx). (8) 

The state propagating between the insertions of the x operators at t\ and t% cannot be 
the ground state, since x switches parity. If t% — t\ is large enough then the first excited 
state will dominate, by the same argument as used above for ground state domination. 
So, if we can evaluate the ratio of path integrals, we can determine the energy splitting 
E\ — Eq between the first excited and the ground state. The evaluation is very simply 
done by taking an ensemble of configurations generated by the Metropolis algorithm 
above and 'measuring' on each one the value of x(t2)x{t\). The ensemble average of 
this quantity, denoted << x{t 2 )x[ti) >>, is then the ratio above. 

The evaluation will suffer from a statistical error depending on how many config- 
urations are in the ensemble. This will improve as the square root of the number of 
configurations provided that the results on each configuration are statistically indepen- 
dent. Since the configurations were made in a sequence, this will not be strictly true 
and results on neighbouring configurations in the sequence will be correlated. There are 
various statistical techniques, such as binning (averaging) neighbouring results and then 
recalculating the statistical error, that will uncover correlations. If the statistical error 
grows with the bin size then the results are correlated and the binned results should be 
used rather than the raw results. It might also be true that results on some number of the 
initial configurations of an ensemble have to be discarded because these configurations, 
and measurements on them, are not yet typical of the e~ s distribution. For this example, 
the results can obviously be improved statistically by moving t\ and t% along the time 
axis, keeping t 2 — t\ fixed, and averaging those results. 

To extract E\ — Eq we should fit the results as a function of t 2 — 1\ to the exponential 
form of Equation[8] Since t 2 — ti will not be infinite, we can take account of contam- 
ination from higher states by fitting to a sum of exponentials in which we constrain the 
higher E n > E\ . Finally we may also improve the accuracy on the energy by modifying 
the operator x(ti)x(t 2 ) to a product of functions of x(ti y2 ) which has optimum overlap 
with the first excited state and minimum overlap with higher excited states. This will 
allow Equation[8]to become true for smaller values of t 2 — t\ and improve the extraction 
of E\ — Eq from the fits. Indeed we can fit simultaneously to the results from several 
different operators and this again will improve the accuracy with which Ei — Eq can be 
determined. It is a useful exercise to do the calculation above for, say, the simple har- 
monic oscillator potential, V = x 2 jl and m = 1. All the pitfalls of 'measurement' and 
the improvements above have their mirror in lattice QCD calculations, as we shall see, 
but this simple example demonstrates them very clearly in a setting where computer time 
is not an issue. 

An important issue is that of the systematic errors, called discretisation errors, in- 
troduced because of the non-zero lattice spacing. We can make these explicit by rewrit- 
ing the finite difference in terms of the continuum (ie continuous space-time of the real 
world) derivative. The exponential of the continuum derivative is the translation operator 
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for moving from one site to the next. So 

x l+ i -Xi e 9a - 1 ad 2 

= = d+— . (9) 

a a I 

This shows that this form of the finite difference has errors linear in a. The actual size of 
the errors at a given value of a would depend on the effective size of d 2 in the quantity 
being calculated. To halve the errors requires half the lattice spacing at (at least) double 
the computer cost. 

The situation is significantly improved by using an improved discretisation. Axi — 
(xi + i — Xi-%)/2a has errors first at 0(a 2 ). Further improvement is obtained by correct- 
ing for the a 2 error using 

Axi - ^(A) 3 ^ (10) 
6 

where (A) 3 is a discretisation of the third order derivative. Errors are then 0(a A ) and this 
means very much smaller errors at a given value of a, or a huge reduction in computer 
cost for a given error by being able to work at a larger value of a. The improved dis- 
cretisation costs a little in computer power to implement but this is completely negligible 
compared to the computer cost saving of the improvement and it is this that has made 
lattice QCD calculations tractable. 



2.2 Gluon fields in lattice QCD 

Now we have all the ingredients for lattice QCD. The position operator as a function of 
time is replaced by the quark and gluon field operators as a function of four-dimensional 
space-time. We discretise the space-time into a lattice of points (Figure^ in order to be 
able to calculate the Feynman Path Integral numerically, using Monte Carlo methods on 
a computer. The ground-state which dominates the path integral is the QCD vacuum and 
we will be interested in excitations of this which correspond to hadrons. 




Figure 2. A 2-dimensional rendition of a 3 -dimensional cubic lattice. Lattice QCD 
calculations use a 4-dimensional grid. 
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It is important to remember that lattice QCD is not complicated. It is a straightforward 
simulation of the theory on a computer, but a numerically intensive one. It is the theory 
of QCD that is being simulated, not a model. However, we are limited in the things 
that we can calculate. There are also statistical and systematic errors associated with the 
calculations and it is important to understand where they come from, so that you can 
assess the usefulness of a particular calculation for your needs. 

Consider the operator O = (tpip) y (ipip) x . This creates a hadron at a point x and 
destroys it at a point y. This is the QCD generalisation of the x(t\)x(t2) operator of the 
previous section. Then the matrix element in the vacuum of the operator is given in path 
integral form as: 



The path integral runs over all values of the quark and gluon fields ip and A at ev- 
ery point in space-time. Discretisation of space-time onto a lattice makes the number of 
space-time points (and therefore field variables) finite. Continuous space-time (x, t) be- 
comes a grid of labelled points, (xi,ti) or (n^a, n t a) where a is the lattice spacing (this 
doesn't have to be the same in all directions but usually is). The fields are then associated 
only with the sites, ip(x,t) — > ip(n i} n t ). The action must also be discretised, but, as in 
the previous section, this is straightforward, replacing fields by fields at the lattice sites 
and the derivatives by finite differences of these fields. The integral over space-time of 
the Lagrangian becomes a sum over all lattice sites: ( J d 4 x — > J2 n fl4 )' anc ' tne P at ^ m " 
tegral becomes a product of integrals over each of the fields, to be done by Monte Carlo 
averaging. These integrals will be finite, unlike such integrals in the continuum, because 
the lattice provides a regularisation of the theory. The lattice spacing provides an ultra- 
violet cut-off in momentum space since no momenta larger than it /a make sense (since 
the wavelength is then smaller than a). 

To discretise gauge theories such as QCD onto a lattice in fact requires a little ad- 
ditional thought to that described above because of the paramount importance of local 
gauge invariance. The role of the gluon (gauge) field in QCD is to transport colour from 
one place to another so that we can rotate our colour basis locally. It should then seem 
natural for the gluon fields to 'live' on the links connecting lattice points, if the quark 
fields 'live' on the sites. 

The gluon field is also expressed somewhat differently on the lattice to the continuum. 
The continuum is an 8 -dimensional vector, understood as a product of coefficients A b ^ 
times the 8 matrices, Tj,, which are generators of the SU(3) gauge group for QCD. On the 
lattice it is more useful to take the gluon field on each link to be a member of the gauge 
group itself ie a special (determinant =1) unitary 3x3 matrix. The lattice gluon field is 
denoted U^(ni,n t ), where [i denotes the direction of the link, n^, n t refer to the lattice 
point at the beginning of the link, and the colour indices are suppressed. We will often 
just revert to continuum notation for space-time, as in U^x). The lattice and continuum 
fields are then related exponentially, 




-s 



(11) 



(12) 
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where the a in the exponent makes it dimensionless, and we include the coupling, g, for 
convenience. If (x) is the gluon field connecting the points x and x+ 1 M (see Figure[3}, 
then the gluon field connecting these same points but in the downwards direction must 
be the inverse of this matrix, U~ (x). Since the U fields are unitary matrices, satisfying 
WU = 1, this is thenUtix). 



x+1 



x+l 
O 



U-^x + l) = U„ 



UUx) 



Figure 3. The gluon field on the lattice. 

This form for the gluon field makes it possible to maintain exact local gauge invari- 
ance on a lattice. To apply a gauge transformation to a set of gluon fields we must specify 
an SU(3) gauge transformation matrix at each point. Call this G(x). Then the gluon field 
Ufj,(x) simply gauge transforms by the (matrix) multiplication of the appropriate G at 
both ends of its link. The quark field (a 3-dimensional colour vector) transforms by 
multiplication by G at its site. 



Uj?\x) 

v> (9) (V) 

ip (9 \x) 



G{x)U^{x)G\x + l^) 

G(x)i/j(x) 

^{x)tf(x). 



(13) 



To understand how this relates to continuum gauge transformations try the exercise of 
setting G{x) to a simple U(l) transformation, e ~ la ( x )^ and show that Equation 1131 is 
equivalent to the QED-like gauge transformation in the continuum, A 9 = An — d^a. 



X2 




Figure 4. A string of gluon fields connecting quark and antiquark fields (left) and a 
closed loop of gluon fields (right). 

Gauge-invariant objects can easily be made on the lattice (see FigurefU out of closed 
loops of gluon fields or strings of gluon fields with a quark field at one end and an 
antiquark field at the other, eg tp(xi)U fl (xi)U v (xi + 1 M ) . . . U 6 (x 2 — le)^^)- Under 
a gauge transformation the G matrix at the beginning of one link 'eats' the G* at the end 
of the previous link, since G^G = 1. The G matrices at x\ and x 2 are 'eaten' by those 
transforming the quark and anti-quark fields, if we sum over quark and antiquark colours. 
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The same thing happens for any closed loop of Us, provided that we take a trace over 
colour indices. Then the G at the beginning of the loop and the G^ at the end of the loop, 
the same point for a closed loop, can 'eat' each other. 

The purely gluonic piece of the continuum QCD action is 

Scont = / d^X^TlF^F^ (14) 

where is the field strength tensor, 

F MV = dfj,A u - d v Ap + ig[A^ A v ]. (15) 
The simplest lattice discretisation of this is the so-called Wilson plaquette action: 

S la tt=/3^(l-^Re{Tr;y p }); f3 = p. (16) 

In fact the 1 here is irrelevant in the lattice QCD calculation, giving only an overall 
normalisation that vanishes from the ratio of path integrals, and so it is often dropped 
from Siatt- 



* J 

x 

Figure 5. A plaquette on the lattice. 

Up is the closed lxl loop called the plaquette, an SU(3) matrix formed by multiply- 
ing 4 gluon links together in a sequence. For the plaquette with corner x in the i,j plane 
we have (Figure 

U p ^(x) = Ui{x)Uj{x + U)Uj(x + lj)u]{x) (17) 

Tr in Si at t denotes taking the trace of U p ie the sum of the 3 diagonal elements. S\ a tt sums 
over all plaquettes of all orientations on the lattice. (3 is a more convenient version for 
the lattice of the QCD bare coupling constant, g 2 . This is the single input parameter for 
a QCD calculation (whether on the lattice or not) involving only gluon fields. Notice that 
the lattice spacing is not explicit anywhere, and we do not know its value until after the 
calculation. (This is a difference from the quantum mechanical example of the previous 
section where we had to choose and input a value for a.) The value of the lattice spacing 
depends on the bare coupling constant. Typical values of (3 for current lattice calculations 
using the Wilson plaquette action are (3 w 6. This corresponds to a w O.lfm. Smaller 
values of (3 give coarser lattices, larger ones, finer lattices. This is obvious from the 
asymptotic freedom of QCD, which tells us that the coupling constant g 2 goes to zero 
at small distances or, equivalently, high energies. (3 is the inverse of the bare coupling 
constant at the scale of the lattice spacing and therefore tends to oo as the lattice spacing 
goes to zero. 



Lattice QCD 



9 



That Siatt of Equation[n[]is a discretisation of Scant is not obvious, and we will not 
demonstrate it here. It should be clear, however, from Equations ^] and ^] that Si at t 
does contain terms of the form d^A^ needed for the field strength tensor, F^ v . 

SWt is gauge-invariant, as will be clear from earlier. Thus lattice QCD calculations 
do not require gauge fixing or any discussion of different gauges or ghost terms, as would 
be required for continuum calculations using perturbation theory. Our lattice calculation 
is fully non-perturbative since the Feynman Path Integral includes any number of QCD 
interactions. In contrast to the real world, however, the calculations are done with a 
non-zero value of the lattice spacing and a non-infinite volume. In principle we must 
take a — > and V — > oo by extrapolation. In practice it suffices to demonstrate, with 
calculations at several values of a and V, that the a and V dependence of our results is 
small, and understood, and include a systematic error for this in our result. 

Before discussing quarks we illustrate here how a lattice calculation is done with only 
gluon fields, how the lattice spacing is determined, and what the systematic errors are. 

A quantity that can be calculated in the pure gluon theory is the expectation value of 
a closed loop of gluon fields. In fact this can be related to the QCD potential between 
an infinitely massive quark and antiquark. Although not directly a physical quantity, this 
is something about which we have some physical understanding. An infinitely massive 
quark does not move in spatial directions and so simply generates a path which is a string 
of gluon fields in the time direction. If this is joined to a string generated by an antiquark 
a distance R in lattice units away then a rectangular RxT loop of gluon fields is created. 
This is known as a Wilson loop, see Figure[4]right. 

To 'measure' this in pure gluon QCD, we generate configurations of gluon fields 
with probability e~ Slatt where Si a tt is a discretisation of the pure gluon QCD action 
given by Equation^] On each of these configurations we calculate the RxT Wilson 
loop, averaging over all positions of it on the lattice. We then average over the results on 
each configuration in the ensemble to obtain a final answer with statistical error, for lots 
of values of R and T, We have 



Typically we need many hundreds of configurations in an ensemble for a small statistical 
error at large R and T. 

The ensemble average of O is related to the heavy quark potential by similar argu- 
ments to those used for the operator x(t\)x(ti) in section 2.1. One end of the Wilson 
loop creates a set of eigenstates of the Hamiltonian that are based on a massive quark- 
antiquark pair. These eigenstates are produced with different amplitudes by the Wilson 
loop operator and have different energies. In this case there is no kinetic energy, so the 
energies are those of the heavy quark potential. The different eigenstates propagate for 
time T and, if T is large, the ground state eventually dominates. 



By fitting the results as a function of the time length, T, in lattice units, the heavy quark 
potential in lattice units, aV(R), is obtained. V' is some kind of excitation of the poten- 
tial which we will not be interested in here. The heavy quark potential at short distances 




«0»=Ce~ 



*V{R)T + c , e 



~aV'(R)T 



+ 



(19) 
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should behave perturbatively and take a Coulomb form. At large distances we expect a 
'string' to develop which confines the quark and antiquark and gives a potential which 
rises linearly with separation. We can therefore fit the lattice potential to the form 

4 n (r\ 

aV{r = Ra) = ^ + aa 2 R + C (20) 

3 R 

where C is a 'self-energy' constant that appears in the lattice calculation. If the results 
for aV(R) are plotted against R, the slope at large R is the 'string tension', a, in lattice 
units, ie aa 2 . Phenomenological models of the heavy quark potential give values for y/a 
of around 440 MeV. Using this value for a and the result from the lattice of aa 2 , gives 
a value for a. This is often quoted as a value for a -1 in GeV. Note that a -1 in GeV 
= 0.197/(ainfm). 

6 
4 
2 

O 

K 

g 
> 

-2 
-4 
-6 

12 3 4 

R/R n 

Figure 6. The heavy quark potential in units of the parameter, ro, as a function of 
distance, r, also in units o/ro. The calculations were done in quenched (pure glue) QCD 
at a variety of different values of the lattice spacing, corresponding to the different values 
of (3 quoted. (Bali 2000) 

The results for aV can then be multiplied by a -1 to convert them to physical units of 
GeV and, having removed the constant C, the results can be plotted as a function of the 
physical distance, r in fm. If discretisation errors are small, then results at different values 
of the lattice spacing should be the same. Figure [5] shows results for the heavy quark 
potential on relatively fine lattices at different values of j3 using the Wilson plaquette 
action for S g (Bali 2000). V(r) is not in fact given in GeV here, nor is r in fm, but both 
are given in terms of a parameter called ro (Sommer 1994). This is the value of r at 
which r 2 dV/dr = 1.65 and is a commonly used quantity to determine the lattice spacing 
(or at least relative lattice spacings), rather than the string tension, ro is not a physical 



X = 6.0 
: + (8=6.2 
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parameter, and as such is not available in the Particle Data Tables. We shall see later that 
there are good hadron masses to use for the determination of a and these can also be used 
to determine the value of tq. Meanwhile, ro ~ 0.5fm. The results at different values of 
P lie on top of each other and this gives us confidence that the discretisation errors are 
small. 

• • • • f f f f f f f 



• • • • 



• • • • 



• • • • • • • • • • • 

Figure 7. A physical object on two lattices of different lattice spacing. On the right our 
updating algorithm takes much longer to register a significant change on the length scale 
of the object. 

The values of (3 used correspond to rather fine lattices. The coarsest lattice is at j3 
= 6.0 and has a=0.1fm. The finest lattice is at p = 6.8 and has a = 0.03fm. The finest 
configurations are very expensive to generate. Yet if we looked closely at the results at 
P = 6.0 we would be able to see discretisation errors at the few percent level. If we want 
accurate results on lattices that are coarse enough for affordable calculation (especially, 
as we shall see, once we include quarks) then we must improve the discretisation of 
the action. The cost of lattice calculations grows naively as a~ 4 because of the four- 
dimensional space-time. In fact it is even worse than this because of critical slowing- 
down. Physical distances on the lattice grow in lattice units as the lattice spacing gets 
smaller (see Figure^}. This means that algorithms that update configurations by making 
local changes to the fields, get slower and slower at making a change on a distance scale 
of r as a is reduced. Then the cost of the calculation of, say, V(r) actually grows as 
a~ 6 . It becomes imperative to improve the action, rather than to attempt to beat down 
the discretisation errors by reducing a. 

Improving the gluon action is at first sight straightforward. We expand the Wilson 
plaquette action in powers of a and notice that it has a 2 errors when compared to the 
continuum QCD gluon action. We add a higher dimension operator to the action in the 
form of a 2 x 1 Wilson loop to cancel this error: 



Ci and c 2 are chosen to remove the a 2 error and naively would have the value 1 . How- 
ever, in a quantum field theory like QCD, the value of c\ y 2 becomes renormalised by 
radiative corrections. S g .i a tt must reproduce S gtCon t to a required level of accuracy and 
so ci,2 must absorb the effect of the differences in gluon radiation between the continuum 




(21) 
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and lattice versions of QCD. This difference arises from gluon radiation with momentum 
larger than 7r/a which does not exist on the lattice. Gluon radiation at these high mo- 
menta is perturbative and so we can calculate c\ t i as a perturbative power series in a s . 
Provided a is small enough so that a s (ir/a) is small enough, the C1.2 will not be very 
different from 1 . There is a complication, however, and that arises from the way in which 
the lattice gluon field is related to the usual gluon field A^. The exponential rela- 
tionship in Equation^]means that the lattice perturbation theory contains vertices with 
many powers of An- Although these are suppressed by powers of g, they produce rather 
large contributions and have to be taken into account. In fact they appear as so-called 
'tadpole' diagrams which take the same form in many different processes (Lepage and 
Mackenzie 1993). This allows them to be substantially removed in a universal way by 
the simple expedient of estimating how far the gluon link field, U^, is from the value 1 
(the 3x3 unit matrix) that it would take in the continuum, g 2 — ► limit. We can measure 
this, for example, from the average plaquette (traced and divided by 3 so that it would be 
1 in the continuum limit). This contains 4 U fields so we take the fourth root to determine 
the 'tadpole-factor', uq. Taking 



is called 'tadpole-improvement'. If this is done then the calculation of C12 does indeed 
give well-behaved perturbative expansions which do not differ substantially from the 
naive value of 1 . The improved gluon action becomes 



dropping the constant term in the action. Taking ci,2 to be 1 is the tree-level tadpole- 
improved Symanzik action and this gives results significantly better than the Wilson pla- 
quette action on coarse lattices (Alford et al. 1995, Lepage 1996). Note that (3 here 
cannot be compared directly to that for the Wilson plaquette action. It only makes sense 
to compare results between different actions at the same value of the lattice spacing. 

Figure [8] compares results for the heavy quark potential, calculated as discussed 
above, for the two actions on coarse lattices with lattice spacing 0.4fm (Alford et al. 
1995). The results with the Wilson plaquette action show clear discretisation errors in 
the fact that V(r) is not a smooth curve, ie it does not have rotational invariance. It re- 
flects the fact that distances on the lattice are not in general travelled in straight lines and 
the result will depend on the actual path used between two points. The curve obtained 
with the improved gluon action is much smoother. As discretisation errors are removed, 
the path dependence of the distance measurement is reduced and the rotational invariance 
of the continuum is restored. Indeed C12 could be fixed non-perturbatively (ie within the 
lattice simulation) by demanding rotational invariance and it is clear from this that results 
close to 1 after tadpole-improvement would be obtained. 

In section 3 we will describe recent lattice results and these use an improved gluon 
action to reduce the discretisation errors. In fact the action is improved beyond that used 
in Figure [8]by including the first a s terms in c\ and C2 and an additional operator that 
appears at 0(a s ) made of a 3-dimensional parallelogram of U fields (Alford et al. 1995, 
Lepage 1996). 
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Figure 8. The heavy quark potential calculated in pure gluon QCD at a lattice spacing, 
a, ofOAfm. The left plot uses the Wilson plaquette action, the right plot the tree-level 
tadpole-improved Symanzik action. The right plot has a much smoother curve than the 
upper one reflecting the improvement of discretisation errors (Alford et al. 1995). 



2.3 Quark fields in Lattice QCD 

Quarks represent a big headache for Lattice QCD. Because quark fields anticommute, 
they cannot be handled using ordinary numbers on a computer. We therefore have to 
perform the integral over the quark fields in the path integral by hand. In fact this is easy 
to do because of the way that the quark fields appear in the QCD action. Including the 
Dirac piece of the action for quarks, we have 

Z = j VUVtpV^e- SQCD 
Sqcd = S g + J2^-D + m)^ = S g +^M^. (24) 

X 

The quark field, ip is a 12-component (3 for colour and 4 for spin) field at every lattice 
point, so M is a matrix with 12 x L 3 T rows and columns for an L 3 T lattice. D is a 
covariant derivative that includes a coupling with the gluon field. The integral over quark 
fields gives 

Z = j VUdet{M)e- s o- QCD = J VU e~ §QCD (25) 

where S = S 9 ^qc d + In det(M), one detM factor per quark flavour. 

Typically in lattice QCD we want to calculate the mass of a hadron made of quarks 
and antiquarks. To do this we must calculate the expectation value of a product of appro- 
priate operators made of quarks and antiquarks. A suitable operator that creates a meson 

is ^ a ' Q ' Jl r Q! ' 3 V a, ' 3 " f2 where Y is some combination of 4 x 4 7 matrices that give the 
right spin-parity ( J p ) quantum numbers to the meson, a is a colour index, a and (3 are 
spin indices and fa are flavour indices. We include the conjugate operator to destroy the 
meson T lattice spacings away in time. Then, suppressing colour and spin, 
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On integration, the right-hand-side becomes (for the case where /i ^ /2): 

/ VU Ti- ap i„,coiour,jf(M-| 0iT) rM^ J T0) r)detM e - s ^ CD 

J X>f7 detM e- s 9.Q^ ' (27) 

This calculation requires us to generate sets of gluon field configurations with probability 
e~ SQCD , calculate the trace over spin and colour of the M~ x factors on each configura- 
tion and then average these over the ensemble. The calculation is illustrated in Figure|5] 
It is known as a 2-point function calculations because there are two operators at times 
and T, represented by the filled ovals. The straight lines at the top and bottom indi- 
cate the valence quark propagators (the M _1 factors above that connect the creation and 
annihilation operators for that particular valence quark). As the quark propagates it in- 
teracts any number of times with the other quark through the lattice gluon field and this 
is indicated by the curly lines. Some of these gluon lines may include the effect of a sea 
quark-antiquark pair (as on the left of the diagram). 




Figure 9. An illustration of the calculation of a simple 2-point function for a hadron in 
lattice QCD. 

If we have f\ = f% for a flavour-singlet quantity then there are additional 'discon- 
nected' pieces where the M^ 1 factors are generated from each operator separately. This 
gives, for example, Tr(M^~ 1 j)Tr(M^, 1 r - ) ). These are very difficult to calculate, being 
very noisy, and accurate results for light flavour-singlet mesons are not yet available. For 
heavy mesons, cc and bb, this is expected to be a very small effect and is usually ignored. 

The hadron mass is determined from the results through the usual multi-exponential 
form 

^(\0H\T)H(0)\0) =«Tr spin , coloul , if (M / - 1 rAf / r 1 r) »=Y J C n e- m - aT . (28) 

n 

n runs over radial excitations of the hadron with a particular set of J p and flavour quan- 
tum numbers, and m n a are the different hadron masses. In fact the formula above is 
only correct for a lattice of infinite time extent in lattice units, T'. On the finite lattices 
we must use there are additional terms from the possibility of quarks 'going backwards' 
round the lattices. These terms take the form e -m ™ a ( T -T ) and can readily be taken ac- 
count of in the fit, although their form does depend on the hadron being studied. C n in 
the above equation is related to the square of the matrix element < 0|if |n > by analogy 
with Equation[8]in section 2.1. The size of C n then depends on the form of the operator 
H used to create and destroy the hadron. Any operator with the right quantum numbers 
can be used and if we make a set of such operators we can calculate a whole matrix of 
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correlators and fit them simultaneously for improved precision. Typically we use opera- 
tors made of quark and antiquark fields but separated in space according to some kind of 
'wavefunction', (known as a 'smearing') eg il) x <j>{x — y)4>y To make a gauge-invariant 
operator of this kind either requires strings of Us to be inserted between x and y or for 
the gluon field configuration to be gauge-fixed, typically to Coulomb gauge. It may be 
possible to choose <j> so that a particularly good overlap with one of the states of the sys- 
tem (eg the ground state) is obtained, and poor overlap with the others. Then Co would 
be large in the equation above, and the other C n small, and this would give improved 
precision for m a. The masses of the radial excitations are of interest too and smear- 
ings which improve overlap with them have also been studied. With particular operators, 
H, C n contains information of physical relevance. For example the matrix element of 
the temporal axial current, Ja between the vacuum and a pseudoscalar meson gives its 
decay constant, related to the rate for its decay purely to leptons via a W boson. The 
calculation in lattice QCD of matrix elements of this kind is discussed in section 3.3. 

The way in which a lattice QCD calculation with quarks is done makes a clear dis- 
tinction between valence quarks and sea quarks. The valence quarks are those that give 
the hadron its quantum numbers and these give rise to the A/ -1 factors in Eauationl27l 
The sea quarks are those that are produced as quark-antiquark pairs from energy fluctu- 
ations in the vacuum (quark vacuum polarisation). They give rise to the detA/ factors in 
Equation^] The important sea quarks are the light ones, u, d and s. The heavy quarks 
c and b have no effect as sea quarks and the t quark does not even have to be considered 
as a valence quark since it does not form bound states before decaying. 

Manipulations of the matrix M are computationally costly. Even though it is a sparse 
matrix (with only a few non-zero entries) it is very large. There are various computational 
techniques for calculating the M _1 factors and the detM factor is included by repeated 
determination of the calculation of M _1 . This makes the inclusion of detAf very costly 
indeed. It becomes increasingly hard as the quark mass becomes smaller because M 
becomes ill-conditioned. (The eigenvalues of M range between some fixed upper limit 
and the quark mass, so this range increases as m q — > 0). In the real world the u and 
d quarks have very small mass and so in the past there has not been sufficient computer 
power available to include them as sea quarks or, if they have been included, their masses 
have been much heavier than their real values. 

Missing out sea quarks entirely is known as the quenched approximation. It is 
clearly wrong, but for a long time the presence of other systematic errors and poor statis- 
tics obscured this fact. More recently it has become clear that the systematic error in 
the quenched approximation is around 10-20%. When light quark vacuum polarisation 
(detA/) is included the calculation is said to be 'unquenched' or 'dynamical'. The sea 
quarks are then also called dynamical quarks. We will see in the results section that it 
is now possible to include realistic quark vacuum polarisation effects and the quenched 
approximation can be laid aside at last. 

Our earlier discussion makes clear that the computational cost of including light 
quark vacuum polarisation will require a very good (ie highly improved) discretisation 
of both the gluon and quark actions. We have discussed the gluon action earlier. For the 
quark action there are several possibilities, or formulations. We will concentrate here on 
those formulations that have already been used in unquenched simulations. 

The naive quark action is a straightforward discretisation of the Dirac action in Eu- 
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clidean space-time: 



S q = y~] ^(x)(7 • A + ma)ip(x). 



(29) 



X 



The finite difference A is 



(30) 



showing explicitly how the gluon fields appear, coupled to the quarks, to maintain gauge 
invariance in the action. QCD with quarks has parameters, in addition to the coupling 
constant, that are the quark masses. In lattice QCD these appear as quark masses in lattice 
units, ma. As with the lattice spacing/coupling constant, the quark masses are not known 
a priori and must be adjusted as a result of a lattice calculation. We therefore iterate until, 
for example, we get a particular hadron mass correct. Since the quark mass is in lattice 
units, finer lattices will require smaller values of ma than coarser ones. 

The naive quark action has several good properties. It has discretisation errors at 
0(a 2 ). The finite difference piece is anti-hermitian and therefore has purely imaginary 
eigenvalues that appear as ±iA. With the mass term the eigenvalues of M are then 
ma ± i\. This is the same as in the continuum and follows from the chiral symmetry of 
the action ie our ability to rotate independently left- and right- handed projections of the 
quark field. It means that the bare quark mass is simply multiplicatively renormalised 
and the renormalised quark mass goes to zero at the same point as the bare quark mass. 
This is important because we need to work with very small (renormalised) quark masses 
for the u and d quarks and we want to be able to go in that direction by just taking ma to 
be smaller and smaller (if we have enough computing power). 

Chiral symmetry is spontaneously broken in the real world giving rise to a Goldstone 
boson, the ir, whose mass consequently vanishes at zero quark mass (m 2 oc m q ). For 
small quark masses, where m^ is small, we have a well-developed chiral perturbation 
theory which tells us how hadron masses and properties should depend on the u/d quark 
masses (or equivalently mj) and we can make use of this to extrapolate down to physical 
u/d quark masses from the results of our lattice QCD simulation provided that we are 
able to work at small enough u/d quark masses to be in the regime where chiral perturba- 
tion theory works. In general m u / d < m s /2 is necessary for an accurate extrapolation. 
We will discuss this further in section 3. (Arndt 2004) 

It would seem that we have everything ready to do the simulations but we must first 
discuss the infamous doubling problem. The naive quark action in fact describes 16 
quarks rather than 1. This is demonstrated most easily in 1 -dimension in the absence 
of gluon fields. It suffices to compare the continuum derivative in momentum space, ie 
p, with the Fourier transform of the finite difference on the lattice, which is s'm(pa)/a. 
These are plotted in Figure ^3 as a function of p between —ir/a and n/a, the limits 
over which p makes sense on the lattice (the first Brillioun zone). Momenta larger than 
7r /a reappear as negative momenta larger than —ir/a so these points are periodically 
connected. Around p w the sine function goes through zero and mimics a straight 
line up to a 2 errors as expected. The problem is that this is also true around p w ir/a 
(with opposite slope). This means that there is a continuuum-like solution of the Dirac 
equation around p — it /a as well as around p = ie there are 2 quarks instead of 1 . In 
4-d (and including 7 matrix algebra as well) the picture is more complicated but the basic 
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Figure 10. (Left) The continuum derivative and lattice finite difference compared as a 
function of momentum within the first Brillouin zone on the lattice. The points p = Ti/a 
and —it /a are joined by lattice periodicity. (Right) Taste-changing interactions for naive 
quarks. A high-momentum (p — tt /a) gluon exchange can change one taste of quark into 
another. 

facts are the same - there are 2 4 quarks instead of 1. The 'doublers' or extra 'tastes' of 
quark live at the edges of the Brillouin zone, where any component of p is close to ir/a. 

This was originally thought to be disastrous and huge efforts, which continue today, 
were made to solve the problem. Wilson introduced a Wilson term into the naive action 
which has the effect of giving mass to the doublers. The Wilson action is then 

x 

4 

= ^U^x)^^ -2^ x + C^(x-l M )^ a _ lM . (31) 

In the absence of gluon fields the Fourier transform of the extra term is 2r sin 2 (pa/ 2) 
which has a maximum at ±ir/a. When folded in to the inverse quark propagator, taking 
account of the 7 matrix structure, this then prevents the inverse propagator from vanish- 
ing at the edges of the Brillouin zone, effectively giving the doublers a large mass which 
increases as a — > so that they are completely removed in the continuum limit. 

This solution of the doubling problem was used for many years. It does have two big 
disadvantages. One is that discretisation errors now appear at 0(a) which causes large 
systematic errors even on relatively fine lattices. The other is that the Wilson term does 
not respect chiral symmetry and so we lose the simple connection between the bare quark 
mass parameter and the renormalised quark mass and have to search in ma for the point 
where m^ vanishes. The eigenvalue spectrum is now complicated as well which creates 
numerical difficulties. Discretisation errors are ameliorated in the improved form of the 
action known as the 'clover' action. The 0(a) errors are cancelled by a a^F^ term 
which is discretised as a set of plaquettes around a central point that looks like a 4-leaf 
clover. The coefficient of the clover term can be set using tadpole-improvement or non- 
perturbatively. Again it is close to 1 after the tadpole-improvement with a well-behaved 
perturbative expansion. Work using the clover action continues and it is being used 
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for unquenched simulations (Ishikawa et al. 2005). A newer form of the action, called 
twisted mass QCD, is also being developed. This overcomes several of the numerical 
problems and is a promising development for the future (Frezzotti 2005). 

Another relatively new development is that of a set of actions that maintain chiral 
symmetry on the lattice but also solve the doubling problem. These require an enormous 
increase in computer power because they either need a matrix inversion inside a matrix in- 
version to calculate M _1 (overlap quarks) or a calculation in 5-dimensions (domain wall 
quarks). Calculations using these formalisms have then been restricted to the quenched 
approximation because of the cost. In future they may be possible to implement includ- 
ing realistic quark vacuum polarisation effects on very powerful supercomputers (Chiu 
2004). 

The recent progress in lattice QCD calculations has come, however, from returning to 
the naive discretisation of the quark action. In fact the doubling problem is not a problem 
if the doublers are simply copies of each other, because we can then simply 'divide by 
16' inside our calculation. If M has a 16-fold degeneracy in its eigenvalue spectrum then 
taking (detM) 1 / 16 gives us the effect of 1 taste for each flavour. The problem is that the 
different tastes of quarks do interact with each other and change taste and this splits some 
of the degeneracy. The process by which this happens is a high momentum (ir/a) gluon 
exchange as in Figure ^3 which moves the quark from one Brillouin zone boundary to 
another. These taste-changing interactions are an artefact of using the lattice and they 
induce large discretisation errors, even though they are formally 0(a 2 ) (Lepage 1998b). 
To improve the naive action then requires removing the unphysical taste-changing in- 
teractions at leading order as well as the usual a 2 errors from discretising a derivative 
as a symmetric finite difference, discussed in section 2.1. This markedly improves the 
discretisation errors, but also significantly reduces the taste-changing interactions, and 
makes this action a good one for lattice simulations. 

The naive quark formalism is fast numerically because it can be converted to the 
staggered quark formalism in which there is no quark spin. It is a remarkable feature that 
the quark fields can be rotated so that the naive quark action (or its improved variant) is 
diagonal in spin space. We have only to simulate one spin component on the lattice and 
then, if required, we can reconstruct all quantities in terms of naive quarks. If we take 

rl>(x) -> n(x) X (x); $(x) ^ x(x)tf (x) (32) 

where 

3 

the unimproved naive action becomes 

^>(x)(7 ■ A + ma)^ — x(x)(a(x) ■ A + ma)x(x) 

where a(x) is diagonal in spin space, a^(x) — (— \y x o+---'-c^-i _ Rewriting the action 
in terms of a single-spin (but 3-colour) component field x> we nave me unimproved 
staggered quark action: 

Sf =E^^E a ^^W fal f ~ U U X ~ 1 m)Xx-iJ +max x }. (33) 

X fj, 
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The staggered quark action does not solve the doubling problem, but removes an 
exact 4-fold degeneracy from it. There are then 4 tastes for every quark flavour that we 
include with M instead of 16. If we improve the action so that the a 2 taste-changing 
interactions are small then the tastes will be close to being copies of each other. We will 
expect to see a 4-fold close-to-degeneracy in the eigenvalue spectrum of the improved 
staggered M, which will become closer and closer as a — > 0. This has been demonstrated 
numerically (Follana et al. 2004) and justifies 'dividing by 4' by taking (detM) 1 / 4 in the 
path integral to include one quark flavour. The improvement of the staggered quark action 
to remove the leading-order taste-changing interactions of Figure[TO|is done by replacing 
the gluon fields that appear in the finite difference of Equation[33]with a combination of 
gluon fields that remove the coupling between the quark and a single p = ir/a gluon. 
The combination 'smears' out the gluon field in directions perpendicular to its link by 
including in combination with a single U^(x) paths such as 'the staple', U v (x)U fl (x + 
l v )Ul(x + lfj). The simplest improved staggered action with a 2 errors removed uses 
paths up to those containing 7 U fields in all 3 directions perpendicular to a given link 
with tadpole-improvement. It is called the asqtad action and has been successfully used 
in extensive unquenched simulations as described in section 3 on results (Lepage 1998b, 
Orginosef al. 1999). 

All quark formulations must give the same physical results in the continuum a — > 
limit and this can be usefully checked in the quenched approximation. Figure fTTI shows 
results for the mass of a vector light meson (a 'p') at an arbitrary fixed physical quark 
mass as a function of lattice spacing (Davies et al. 2005). Both axes are given in units 
of r±, defined in a similar to way to ro, but r\ w 0.35fm. The x axis is in terms of the 
squared lattice spacing since most of the formalisms plotted have discretisation errors 
at 0(a 2 ) or better. The Wilson formalism, with errors at O(a), is clearly much worse 
than the other formalisms on this plot, with a very steep a dependence, but it does give a 
consistent continuum limit. 

Heavy quarks (b and c) represent a rather different set of issues to those for light 
quarks in lattice QCD. They have large quark masses, tuq, and therefore large values of 
rriQa at any value of a at which we are able to do QCD simulations. This means that we 
risk large discretisation errors when handling these quarks if we use a formalism in which 
the errors are set by the size of ma. If rriQa > 1 then no amount of improvement will give 
a good discretisation for these quarks. This is particularly true for the b quark which has a 
mass of around 5 GeV. To reduce rriQa below 1 requires a -1 > 5 GeV or a lattice spacing 
< 0.04 fm which is incredibly expensive to simulate (as well as being wasteful). Luckily 
physical understanding of these systems comes to our rescue here. The experimental 
spectrum of hadrons containing b and c quarks shows clearly that the masses of these 
hadrons are much larger than the differences in mass between the hadrons in different 
radial and orbital excitations (and these differences are in fact very similar for b and c 
systems). The differences reflect typical kinetic energies and momenta inside the hadrons 
of a few hundred MeV to a GeV, ie of the same size or somewhat larger than those 
typically inside light hadrons but << ttiq. The quark mass itself is not important for 
the dynamics of the bound state but simply provides an overall mass shift (Davies 1998). 
If we formulate the quark action in a way that simulates accurately nonrelativistic quark 
momenta and kinetic energies we will be able to accurately determine the properties of 
the hadrons and it is these scales that will also control the discretisation errors. 
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Figure 11. Results for the vector light meson mass at fixed physical quark mass in 
quenched lattice QCD using a variety of gluon and quark actions (for the valence 
quarks). The mass is given in units of T\ and plotted as a function of a 2 , also in units 
of T\. The lines represent fits to the lattice results, including appropriate discretisation 
errors for each case, that require all lines to have a common continuum limit. A good fit 
is obtained (Davies et al. 2005). 



This means using a nonrelativistic effective field theory and the example we give here 
is Nonrelativistic QCD, NRQCD. The lattice NRQCD action is a discretised nonrelativis- 
tic expansion of the Dirac action, matched to continuum QCD to the desired order in the 
heavy quark velocity, vq and in the strong coupling constant, a s (Lepage et al. 1992). 
The first few terms of the continuum NRQCD quark Lagrangian are: 

C Q =^{D t -^—- Ci ^- + ...)^ (34) 

where we include just the non-relativistic kinetic energy term and the coupling between 
the quark spin and chromomagnetic field. Higher order terms give the spin-orbit interac- 
tion and Darwin term etc. ip is a 2-component spinor here but the antiquark Lagrangian 
is simply related to the quark one and so the antiquark propagator is easily determined 
without a separate calculation. On the lattice the covariant time and space derivatives 
above become finite differences with U fields included. The presence of a single time 
derivative means that the calculation of the quark propagator, M~ 1 , can be solved as an 
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initial value problem on one pass through the lattice. This is much faster than the itera- 
tive methods needed for light quarks. On the lattice tuq becomes mqa, the quark mass 
in lattice units, and this must be determined by getting a heavy hadron mass correct, for 
example the T mass. This is slightly more complicated than in the light quark mass case 
because the direct mass term is missing from the Lagrangian, so the energy exponent of 
an T correlator at zero momentum will not give the mass (although energy differences do 
give mass differences). Instead we have to calculate hadron energies at non-zero spatial 
momentum and extract the mass of the hadron from its kinetic energy. 

NRQCD contains, by design, the low momentum physics of QCD for heavy quarks. 
High-momentum interactions of continuum QCD are missing, and the lattice in any case 
does not allow forp > ir/a. The effect of these missing high-momentum interactions 
is included in NRQCD through a renormalisation of the coefficients, Cj, of the higher 
order terms. Because the effects that we are talking about are high momentum, the 
calculation can be done in perturbation theory as a powers series in a s . Again small 
deviations from 1 are found after tadpole-improvement of the U fields is implemented, 
provided mQa is not large. The Cj contain functions of mga, including inverse powers 
of rriQa, multiplying powers of a s and so we cannot take a to zero or we lose control of 
the convergence of the Cj. This is not a problem in practice since, as discussed earlier, 
discretisation errors can now be controlled in this formalism to high accuracy at values 
of the lattice spacing that we can simulate and there is no need to take a very small. 

For charm quarks the situation is not as clear as for b quarks because m c a w 1 on 
lattices in current use. NRQCD does not necessarily work well and it may be more ad- 
vantageous to take a relativistic formalism of the c quark and attempt to improve it highly 
to reduce the discretisation errors in m c a to a reasonable level. The Fermilab formalism 
does a mixture of both these things by using a relativistic action, the clover action, but in- 
terpreting it non-relativistically where necessary to avoid the largest discretisation errors 
(El-Khadra et al. 1997). This is the method currently being used most extensively for c 
quarks and results will be described in section 3. Another promising method is to use 
anisotropic lattices in which the lattice spacing in the time direction is much finer than 
that in spatial directions. This then allows m c at to be small without requiring a huge 
amount of computer power to make a fine grid in 4-dimensions. The price of this is a 
relatively complicated action, however. 



3 Results 

Recent lattice results which include for the first time realistic quark vacuum polarisation 
effects, have changed the landscape of lattice calculations. We will concentrate on those 
results here and the possibilities for the future that they engender. 

The new results are based on gluon field configurations generated by the MILC col- 
laboration and analysis by the HPQCD and MILC collaborations (Davies et al. 2004). 
They have worked with a highly improved gluon action and with improved staggered 
(asqtad) quarks, as described in the previous section. The effect of the quark vacuum 
polarisation is included in the generation of the gluon configurations by taking the fourth 
root of the determinant of the improved staggered quark matrix for each flavour of quark 
included in the sea. Three flavours of quarks are included in the sea, u, d and s. This is 
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believed to be all that is necessary, since a perturbative analysis shows that the effects of c 
or b quarks in the sea is very small (Nobes 2005). The u and d quarks are taken in fact to 
have the same mass (as in all current lattice calculations) because this provides for some 
simplification and should give only a small error when comparson is made to appropri- 
ately isospin-averaged experimental quantities. The configurations are then referred to 
as '2+1 ' flavour unquenched configurations. 

As described above, it is not possible to perform the lattice calculations at the physical 
u/d quark mass (since we take m u = m d this would be {m u ^ ph y S +m d _ phys ) /2). Instead, 
many different ensembles of configurations are made with different values of m u / d and 
extrapolations to find the physical point must be done. These chiral extrapolations will 
be accurate provided that we are close enough to the physical point and we are able to 
constrain the functional form of the dependence on m u / d with chiral perturbation theory. 

It is convenient, when discussing the u/d quark masses used in the lattice calcula- 
tions, to give them in terms of the strange quark mass, m s . This is because there is no 
problem (with current computers) in doing lattice QCD calculations that include the s 
quark vacuum polarisation. However, it is also clear that, in the real world, there is a 
very visible difference between s quarks and u/d quarks. Indeed, if we consider chiral 
perturbation theory as an expansion in powers of x q — (?nps/(47r/ T ) 2 where mps is 
the mass of the pseudoscalar light meson made from the light quarks, then x s = 0.33 
and x u / d = 0.03. We would therefore not expect low order chiral perturbation theory to 
work well for light quarks with a mass equal to that of the s quark. This is not a problem 
since we do not need chiral perturbation theory to reach the s quark mass. However it 
does mean that our u/d quark mass must be well below m s to expect to be able to reach 
the physical point by chiral extrapolation in m u / d . However, with a set of ensembles 
with m u j d < m s /2 and using next-to-leading order chiral perturbation theory in m u / d , 
it should be possible to perform extrapolations with errors at the few % level. This is 
what has now been done using the configurations generated by the MILC collaboration. 
The MILC configurations range in m u / d down to m u / d = m s /8 which is within a factor 
of three of the real world. This is a huge improvement over previous calculations which 
had only two flavours of sea quarks, meant to represent u and d, but with m u / d > m s /2, 
and it is this improvement that is responsible for the quality of the new results. 

In addition to including the effect of quark vacuum polarisation with light u/d quarks 
and s quarks, the MILC configurations come in sets with three different values of the lat- 
tice spacing. In each set the bare coupling constant must be adjusted as the quark masses 
are changed to keep the lattice spacing approximately the same (as described earlier the 
lattice spacing is determined accurately after the simulation, but preliminary results en- 
able this approximate tuning of the lattice spacing to be done). Having a set of configura- 
tions which include the effect of sea quarks of different light quark mass but which have 
the same lattice spacing is important to allow the chiral extrapolations to be done without 
confusing the (real) effects of changing m u / d with (unphysical) systematic errors from 
having different values of the lattice spacing. The effect of the discretisation errors can 
also be gauged accurately (and an extrapolation to the continuum limit be performed if 
necessary) by having several well-spaced values of the lattice spacing. Again this has not 
been possible with previous lattice calculations. On the MILC configurations the lattice 
spacing values are approximately 0.18 fm (supercoarse), 0.12 fm (coarse) and 0.08fm 
(fine), chosen to be appropriately spaced in a 2 , the expected size of the discretisation er- 
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rors with the improved action used. Most of the results described here will come from the 
coarse and fine sets, since the supercoarse set has only been made recently. A superfine 
set is also now planned, to give a fourth lattice spacing value. In addition another set of 
ensembles is being made with a different bare sea s quark mass. There is no problem in 
principle with working at the correct s mass, but it is hard to tune this accurately until 
after the calculations have been done. On the coarse ensemble the sea s mass is in fact 
25% high and it is 10% high on the fine set of ensembles. Results with a different s mass 
will enable a more accurate interpolation to the correct s mass. 

Whilst singing the praises of the MILC configurations, it should also be said that 
they have a commendably large volume (2.5 fm on a spatial side) and a very long extent 
in the time direction (nearly 8 fm). There are several hundred configurations in each 
ensemble that can be used for calculations. All of these factors lead to small statistical 
errors, indeed smaller than for a lot of calculations that have been done in the quenched 
approximation. The MILC configurations are publicly available at www.nersc.gov. 

3.1 Results on the spectrum 

FigurefT2lsumm arises the new results (Davies et cd. 2004). It shows the ratio of the lattice 
calculation to the experimental number for a range of gold-plated quantities from the tt 
decay constant to radial and orbital excitation energies in the T system, by way of heavy- 
light meson masses and baryon masses. The left-hand panel of the Figure shows results in 
the old quenched approximation and the right-hand panel shows the new results including 
the effect of light quark vacuum polarisation. The clear qualitative difference between 
the two panels is that the left-hand panel shows clear disagreement with experiment for 
a number of quantities because the quenched approximation is wrong. The right-hand 
panel shows instead agreement with experiment for all the quantities within errors of a 
few percent. 

To make Figure [21 the parameters of QCD first had to be fixed. The lattice spac- 
ing was determined from the radial excitation energy in the T system (2S — IS = 
Mr' — My). The u/d quark mass was fixed by finding where the it mass would be 
correct by extrapolation using chiral perturbation theory. The s quark mass was fixed by 
determining where the K meson mass was correct. The c quark mass was fixed using 
the D s and the b quark mass was fixed by getting the T mass correct. All of the hadron 
masses used to determine the parameters are 'gold-plated' ie they have very small decay 
widths and are well below strong decay thresholds. This means that they are well-defined 
experimentally and theoretically and should be accurately calculable in lattice QCD. Us- 
ing them to fix parameters will not then introduce unnecessary additional systematic er- 
rors into lattice results for other quantities. This is an important issue when lattice QCD 
is to be used as a precision calculational tool. 

Having fixed the parameters, we can then focus on other gold-plated masses and 
decay constants and Figure [21 shows the predictions for these other quantities. The fact 
that the right-hand panel demonstrates agreement with experiment for all the quantities 
shown is an indication that the parameters of QCD are unique. Instead of using the hadron 
masses of the previous paragraph to fix the parameters we could have used appropriate 
quantities from Figure[2]and we would have obtained the same answer (and would then 
have been able to predict m^, mx etc). The left-hand panel shows that this is not true in 
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Figure 12. Lattice QCD divided by experiment for a range of 'gold-plated' quantities 
that cover the full range of QCD physics. The unquenched calculations on the right 
show agreement with experiment across the board, whereas the quenched approximation 
on the left gives systematic errors of 0(10%). (Davies et al. 2004) 



the quenched approximation. The results there show disagreement with experiment and 
it is clear, for example, that if we had used the orbital excitation energy in the T system 
(IP — 15) to fix the lattice spacing we would have obtained an answer 10% different. 
The quenched approximation is then internally inconsistent since the parameters depend 
on the hadrons used to fix them. The basic reason is that, as the light quark vacuum 
polarisation is missing, the strong coupling constant runs incorrectly between different 
momentum scales. Therefore hadrons which are sensitive to different momentum scales 
cannot simultaneously agree with experiment. 

In Figure ^]more details are shown of the lattice results for the radial and orbital 
energy levels in the T system in the old quenched approximation and now using the 
MILC configurations with their inclusion of a realistic light quark vacuum polarisation 
(Gray et al. 2003). Our physical understanding of the T system is very good and there 
are a lot of gold-plated states below decay threshold so it is a very valuable system for 
lattice QCD tests and for determining the lattice spacing. We use the standard lattice 
NRQCD effective theory, described above, for the valence b quarks, accurate through v 4 
where v is the velocity of the b quark in its bound state. This means that spin-independent 
splittings, such as radial and orbital excitations, are simulated through next-to-leading- 
order and should be accurate to about 1%. The test of QCD using these splittings is then 
a very accurate one. We do not expect the T system to be very sensitive to the masses 
of the light quarks included in the quark vacuum polarisation, only to their number. The 
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Figure 13. Radial and orbital splittings in the T spectrum from lattice QCD in the 
quenched approximation and including a realistic light quark vacuum polarisation. In 
these plots the b quark mass was fixed from the T mass and the lattice spacing from 
the splitting between the T' and the T. Neither of these masses is predicted. (Top) The 
spectrum of S, P and D levels in the T system obtained from coarse (filled red triangles) 
and fine (open black triangles) quenched lattice calculations and from coarse (filled red 
squares) and fine (open black squares) unquenched calculations. Experimental results 
are shown as lines. (Bottom) Results for different splittings as a function of light u/d 
quark mass. The leftmost points, at lightest u/d quark mass, are the ones included in the 
top plot for the unquenched results. ( Gray et al. 2003) 



momentum transfer inside an T is larger than any of the u, d, or s masses and so we 
expect these splittings simply to 'count' the presence of the light quarks. This lack of 
variation with light quark mass is evident in Figure[T3l 
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Figure 14. Fine structure in the T system. Triangles show quenched results 
(closed red=coarse, open black=fine) and squares show unquenched results (crossed 
green=supercoarse, closed red=coarse, open black=fine). (Top) S-wave splittings be- 
tween T and rjb and between T' and r\' h . The same NRQCD action is used to determine 
the splitting between the B* and the B s and this is shown on the right. There the experi- 
mental situation is summarised by the lines - faint lines show limits for the experimental 
B s results, the darker line shows the more precise result for the B, which is expected to 
be very close to the result for the B s . (Bottom) P-wave splittings between the Po,i 2 Xb 
states, compared to experiment. The results on the right show the lattice prediction for 
the unseen 1 P\ hf, state. (Gray et al. 2003) 
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Figure [3] shows the fine structure in the T spectrum. Since the fine structure is 
determined at leading order by spin-dependent terms that appear first at 0(v A ) this has 
significant systematic errors from missing higher order terms. There are also systematic 
errors arising from missing radiative corrections to the leading v 4 terms. In fact, encour- 
aging agreement with experiment is obtained and we are able to predict the ground-state 
hyperfine splitting between the T and the unseen r/b to be 60(15) MeV. The results for 
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the fine structure can be improved by improving the NRQCD action and this will be 
done in the near future. For comparison to the T splittings, Figure ^] also shows the 
splitting between corresponding mesons made of one b quark and one s quark ie the B s 
heavy-light system. This is calculated on the lattice using the same NRQCD action of 
the valence b quarks and the improved staggered action (as used in the quark vacuum 
polarisation) for the light quark. Once the b quark mass and lattice spacing have been 
fixed from the T system the heavy-light spectrum is entirely predicted by lattice QCD 
and provides another useful test against experiment. 
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Figure 15. Spectrum of mesons containing charm quarks from lattice QCD including 
light quark vacuum polarisation. These results are from the coarse unquenched MILC 
configurations using the T system to determine the lattice spacing and the D s meson 
mass to fix the c quark mass. The Fermilab action was used for the c quark. (Left) The 
charmonium system. (Right) The D s system. (diPierro et al. 2004) 

The charm quark is somewhat too light for good results using the NRQCD action 
above. The results shown in Figure ^] use the Fermilab action described earlier which 
aims to interpolate between light relativistic quarks and the heavy nonrelativistic limit (di 
Pierro et al. 2004). It should be borne in mind that there are fewer gold-plated states in the 
cc spectrum compared to the bb spectrum described above. In fact even the ip' is rather 
close to decay threshold to be sure that its mass is not affected by coupling to virtual 
decay modes which will have a distorted momentum spectrum on a finite lattice volume. 
The right-hand plot of Figure[21shows the D s spectrum. The D s mass was used to fix 
the c quark mass, so that is not predicted. The two newly discovered + and 1 + states are 
shown. The simplest explanation for these states is that they are P-wave states, narrow 
because their masses are too low for the Zweig-allowed decay to D,K. They do decay 
to D s , 7r, however, and are therefore not gold-plated. We expect systematic errors in a 
lattice QCD calculation of their masses, even when light quark vacuum polarisation is 
included. The current lattice results, shown in Figure^] are consistent with this picture 
and certainly do not require any more exotic explanation of the meson internal structure. 

With a good action for b quarks and a good action for c quarks it is possible to predict 
the mass of the B c meson. The results are shown in Figure [16] (Allison et al. 2005) 
compared to the very recent experimental results (Acosta et al. 2005). The quantity that 
is calculated on the lattice is the mass difference between the B c and a combination of 
masses of other mesons containing a b quark and a c quark, enabling some of the sys- 



28 



Christine Davies 



tematic errors to cancel. There are two combinations available on the lattice. The one 
that gives the most accurate result we believe is to calculate the difference between the 
B c mass and that of the average of the T and ip. This is called A^t- It turns out to be 
very small, although non-zero. The difference between the B c and the sum of the masses 
of the D s and B s , Ad b b b , is also useful but less accurate. Both splittings are shown in 
Figure[H)]as a function of the light quark mass included in the quark vacuum polarisation 
on the MILC unquenched configurations. After a study of systematic errors from fixing 
the quark masses and lattice spacing we are able to give a final result for the B c mass of 
6304(20) MeV (Allison et al. 2005). The error is much better than the 100 MeV obtained 
for previous calculations in the quenched approximation. There 100 MeV arose directly 
from the impossibility of consistently determining the parameters of the theory in the 
quenched approximation, a problem that disappears in the new unquenched results. The 
very recent experimental result from CDF is 6287(5) MeV, showing impressive agree- 
ment with the unquenched lattice calculation (Acosta et al. 2005). 




Figure 16. (Left) The differences between the B c mass and, above, the average ofT and 
tj) masses and, below, the B s and D s masses, plotted against the light quark mass in- 
cluded in the quark vacuum polarisation in the MILC unquenched configurations. (Right) 
Predictions for the mass of the B c from a variety of methods. In the centre is the old 
quenched lattice QCD result and on the right the two new lattice QCD predictions. The 
dashed line shows where the B c mass would be equal to the average ofT and %j) masses 
and the broader coloured line is the new experimental result. (Allison et al. 2005, Acosta 
etal. 2005) 

Results for light mesons are equally good on the MILC configurations (Aubin et al. 
2004a, 2004b, Bernard 2001). Figure [TTl shows results for the light meson masses and 
decay constants plotted against light quark mass. The left-hand plot shows the fits to 
chiral perturbation theory for the squared masses of the 7r and K, which must then be 
extrapolated (for the tt) or interpolated (for the K) to find the quark masses which cor- 
respond to u/d (bearing in mind that the u and d quark masses have been taken to be 
the same) and s. Once this has been done, predictions for other light meson quantities 
can be made. Two very clean quantities to calculate are the light meson decay constants, 
/tt or Jk- The decay constant is obtained from the matrix element of the temporal ax- 
ial current between the pseudoscalar light meson and the QCD vacuum and it is related, 
for the charged it and K, to the rate of purely leptonic decay. This rate is well-known 
experimentally for these mesons and can be calculated accurately on the lattice. This is 
because, unlike the matrix elements we will discuss below, the axial current needs no 
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renormalisation for improved staggered quarks. The right-hand plot of Figure [T7l shows 
the chiral extrapolations for the it and K decay constants for the fine MILC lattices. Only 
points for which the u/d quark mass is less than half the s quark mass are used in the fits 
to chiral perturbation theory. The curves, however, show what happens when these fits 
are extrapolated to higher u/d quark masses. For f v it is clear that the heavier u/d quark 
are not on the same chiral perturbation theory curve as the lighter ones ie using heavier 
u/d masses would have distorted the chiral fit and given the wrong result in the chiral 
limit where the u/d quark mass takes the correct value to give the experimental tt meson 
mass. This should be borne in mind when looking at earlier lattice results that have only 
m u / d > m s /2. 




Figure 17. (Left) Chiral fits and extrapolations for light meson ( n and K) meson masses 
from the MILC configurations. The square of the meson mass is plotted against the light 
(u/d) quark mass used in the quark vacuum polarisation. This is denoted mx and given 
in units of the s quark mass used in the quark vacuum polarisation (denoted here m' s ). 
Several different values of valence s quark mass were used so there are several different 
sets of results for the K mass. The lattice results are the points and lines are high-order 
chiral perturbation theory fits. At small light quark mass these are indistinguishable from 
the lowest order straight line. (Right) Chiral fits for f n and Jk, the tt and K decay con- 
stants. The points are lattice calculations on the fine unquenched MILC configurations, 
for two different values of the light quark mass in the quark vacuum polarisation ( fn^ a d ). 
The results are plotted against the valence u/d quark mass in units of the valence s quark 
mass. Note that only points with m u / d < m s /2 (for both valence and sea) were used 
in the chiral fits, shown as lines. The fits were then extrapolated back to compare to the 
lattice results with rn™ l ^ nce > m s /2. (Aubin et al. 2004a, 2004b, Lepage and Davie s 
2004) 

Light baryons are also gold-plated particles and calculating their masses provides a 
further opportunity for predictions from lattice QCD since all the QCD parameters have 
been fixed from the meson sector. Baryons are harder to work with on the lattice and the 
masses are not as precise. The nucleon also requires, for example, a more complicated 
chiral fit and this has not yet been done. Figure IT~8l shows the nucleon results on the 
super-coarse, coarse and fine MILC unquenched configurations (Aubin et al. 2004a). 
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There are signs that the results change slightly with lattice spacing ie there is a visible 
discretisation error. This needs to be studied further. The line shows the low order chiral 
perturbation theory result and the lattice results do seem to be heading towards that on 
the finer lattices. The right-hand plot shows results for the ft baryon made of three 
valence s quarks (Toussaint and Davies 2005). This baryon is not very dependent on the 
u/d quark mass used in the quark vacuum polarisation so chiral extrapolations are much 
simpler. Its mass is sensitive to the s quark mass, however, so the fact that agreement 
with the experimental result is obtained with the s quark mass fixed from the K is another 
strong confirmation that lattice QCD is internally consistent once realistic quark vacuum 
polarisation is included. 

3:0 I 1 1 1 r— n 1 1 — I — 1 — n 3.0 I 1 1 1 1 1 1 1 1 1 1 1 — n 




Figure 18. (Left) The nucleon mass in units of the n parameter from the unquenched 
MILC configurations at three different values of the lattice spacing and for several values 
of the light quark mass, denoted on the x axis by the ratio of it to p masses made of 
the light quarks. The line shows continuum chiral perturbation theory for this quantity 
(Aubin et al. 2004a). The fancy diamond gives the experimental point in these units. 
(Right) The Q baryon mass as a function of sea light quark mass m! divided by sea s 
quark mass. The Q mass is calculated for a valence s quark mass which is the correct 
one, fixed from the K (Toussaint and Davies 2005). The burst gives the experimental 
point. Also on this plot is the splitting between the spin average of the B s and B* and 
the T, a quantity sensitive to the s quark mass, but not the b quark mass. Again it shows 
good agreement with experiment denoted by a burst (Aubin et al. 2004c). 

It is important to realise that accurate lattice QCD results are not going to be obtain- 
able in the near future for every hadronic quantity of interest. What these results show 
is that 'gold-plated' quantities should be now be calculable. Unstable hadrons, or even 
those within 100 MeV or so of Zweig-allowed decay thresholds, have strong coupling 
to their real or virtual decay channel and this is not correctly simulated on the lattice 
volumes currently being used, eg the smallest non-zero momentum on typical current 
lattices exceeds 400 MeV. This could significantly distort the decay channel contribution 
to the hadron mass. Much larger volumes will then be necessary to handle these hadrons. 
Unfortunately the list of non-gold-plated hadrons is a long one - it includes the p, cj>, D*, 
A, N*, glueballs, hybrids etc. Some of these may be more accurately calculable than 
others and qualitative results may also be useful. These points must be borne in mind, 
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however, when making quantitative comparison between lattice QCD and experiment. 



3.2 Determination of the parameters of QCD 

Lattice QCD calculations are an excellent way to determine the parameters of QCD, 
masses and coupling constant, because the method is a very direct one. 
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Figure 19. Quark masses in the M S scheme at a relevant scale (2 GeY for u, d and s 
and their own mass for c and b) as determined from lattice QCD using the unquenched 
MILC configurations (Aubin et al. 2004c, 2004b, Gray et al. 2003, Nobes et al. 2005). 
The lines give the range of the current values quoted in the Particle Data Tables (PDG 
2004). 



For the quark masses we directly determine the bare quark masses in the lattice QCD 
Lagrangian required to give the correct answer for a given hadron mass. With the inclu- 
sion of light quark vacuum polarisation we have seen above that this can be done both 
accurately and unambiguously. Masses are more conventionally quoted in the continuum 
MS renormalisation scheme rather than the lattice scheme. To convert from the lattice 
scheme to MS requires the calculation of a finite renormalisation factor to take into ac- 
count the gluon radiation with momenta larger than n/a that does not exist on the lattice. 
The renormalisation factor can be calculated in perturbation theory since it involves large 
momenta, and it then appears as a power series in a s . This calculation has been done to 
0(a s ) for the light quark masses for the improved staggered quark action (Hein et al. 
2003). The u/d quark mass was fixed from m n and the s quark mass from ttlk as de- 
scribed above. They were then converted to the MS scheme using the renormalisation 
factor. This gave, quoting the masses at the conventional scale: 

m™S{2GeV) = 76(8)MeV, 

mW d (2GeV) = 2.8(3)MeV, (35) 
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where m u /,i is the average of u and d masses (Aubin et al. 2004c). The main source of 
error from the lattice calculation is from unknown higher order terms in the perturbative 
renormalisation factor to convert to MS. A two-loop calculation will be available shortly 
and this will reduce the error to a few %. This is a huge improvement over previous 
determinations of the masses from eg QCD sum rules. Figure[^]shows the results for all 
5 quark masses (u and d are determined separately by matching to different combinations 
of charged and neutral pions and kaons) obtained on the MILC configurations and using 
1-loop matching to convert to the conventional M S scheme, compared to the results from 
the 2004 Particle Data Tables (PDG 2004). It is quite clear that the lattice will take the 
lead in providing accurate quark masses now. 
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Figure 20. (Left) A comparison of the new lattice determination of a s , 0.1170(12) with 
the results from other determinations in the Particle Data Tables (PDG 2004). (Right) 
The determination of a s from the lattice in the V scheme on quenched (rij — 0) and 
unquenched (nj — 3) gluon configurations at various energy scales, d/a for a variety of 
Wilson loop operators. The dashed curves show the expected running of a s in the two 
cases from QCD. (Mason et al. 2005) 

The strong coupling constant, a s , is also well-determined on the lattice (Mason et al. 
2005). The determination of a s proceeds by the calculation in perturbation theory to high 
order of some lattice operator (Trottier 2004). Recently calculations up to and including 
terms in became available for a lot of different Wilson loops and combinations of 
them. These could then be readily 'measured' (ie calculated) on the MILC configura- 
tions. >From each non-perturbative lattice calculation compared to perturbation theory a 
value for a s is extracted at some momentum scale in lattice units. >From the determina- 
tion of the lattice spacing, this scale can then be converted to a physical scale in GeV and 
a s evolved to different scales use the QCD (3 function. Because results at three different 
values of the lattice spacing are available, it is possible to do a consistency check for 
the determination of a s at a given physical scale from three different determinations of 
a particular Wilson loop at three different lattice scales. This allows estimates of fourth 
order terms in the perturbation theory, which puts this calculation into a new regime of 
accuracy for a s determinations. Altogether 28 different loops and loop combinations 
were studied. The a s determined was converted to the A/5 scheme and run to the scale 
of Mz since again this is the conventional comparison point. The final answer obtained 
is 0.1 170(12) which compares very well with other determinations quoted in the Particle 
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Data Tables, see Figurel20l 

Another interesting point is that a s is well able to distinguish between quenched 
and unquenched configurations. Figure |20| shows corresponding determinations of a s 
from MILC quenched and unquenched configurations for different Wilson loops. The 
relevant momentum for the a s determination varies from loop to loop and also depends 
on whether the result comes from the super-coarse, coarse or fine lattices. This enables a 
comparison of a s values with the expected curve for the running. The agreement is very 
good and shows both that the quenched (n f=Q) and unquenched (n f=3) a s figures differ 
markedly and also run differently, both in agreement with the perturbative expectation 
(Mason et al. 2005). 



3.3 Results on matrix elements 

A key point where lattice QCD calculations are needed and can make an impact is in 
the determination of the elements of the Cabibbo-Kobayashi-Maskawa matrix that links 
the quark flavours under weak decay in the Standard Model of particle physics. When a 
quark changes flavour inside a hadron with the emission of a W the quark level process 
is quite simple (see Figure E) , However, this decay unavoidably takes places inside a 
hadron because the quarks are confined by QCD and the QCD corrections to the decay 
rate are significant. This is why the decay matrix element need to be calculated in lat- 
tice QCD so that all of the gluonic radiation around the decay vertex can be taken into 
account. A CKM element Vuf 2 multiplies the vertex but this appears in a simple way in 
the final theoretical answer for the decay rate since there is only one weak vertex in the 
process. A comparison of the experimental decay rate and the lattice QCD results times 

V hh then 8 ivesF /i/2- 

Decay rates which can be accurately calculated in lattice QCD are those for gold- 
plated hadrons in which there is at most one (gold-plated) hadron in the final state. This 
therefore includes leptonic and semi-leptonic decays and the mixing of neutral B and 
K mesons. Luckily there is a gold-plated decay mode available to extract each element 
(except Vtb) of the CKM matrix which mixes quark flavours under the weak interactions 
in the Standard Model: 
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The determination of the CKM elements and tests of the self-consistency of the CKM 
matrix are the current focus for the search for Beyond the Standard Model physics and 
lattice calculations of these decay rates will be a key factor in the precision with which 
this can be done. 
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Figure 21. Current constraints on the vertex of the 'unitarity triangle' made from the 
Cabibbo-Kobayashi-Maskawa matrix. ( CKMfitter 2005) 

Figure|^shows the current 'unitarity triangle' picture in which limits are placed on 
various combinations of elements of the CKM matrix and the result is expressed as a 
search for the vertex of a triangle. The limits that depend on results from B factories 
and the associated lattice calculations, are the dark green ring and the orange and yellow 
rings. The light green hyperbola comes from kaon physics and associated lattice calcula- 
tions. The dark green ring is fixed from semileptonic decays of B mesons to it mesons or 
D mesons. The orange and yellow rings result from mixing of neutral B or B s mesons. 
We will discuss further below the lattice results for these matrix elements. The angles of 
the unitarity triangle, such as sin(2/3), are determined directly by the experiment (light 
blue lines) without theoretical input. 

Decay matrix elements of this kind can be calculated on the lattice from the ampli- 
tudes of the exponentials in the fit functions for hadron masses, Equation[28] For exam- 
ple, the rate at which a B meson decays completely to leptons via a W boson depends on 
the matrix element of the heavy-light axial vector current, = iphJ^sipu, between a 
B meson and the (QCD) vacuum. The matrix element of this current is parameterised by 
the decay constant, fs, so that 

(0\J A JB)= P J B (36) 

(using a relativistic normalisation of the state). For a B meson at rest only Ja (i is relevant 
and the right-hand-side of the equation above becomes to_b/b. 

The way in which Jb is calculated on the lattice is illustrated in Figure |22] On the 
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Figure 22. (Left) The calculation on the lattice of the 2-point correlator for a B meson. 
(Right) The calculation on the lattice of a 2-point function in which a B meson decays 
leptonically. 

left-hand side is an illustration of the usual operator (correlator) whose expectation value 
we calculate to determine the B meson mass or energy. Then, as earlier, 

« H\T)H{Q) »= C n e- E " aT . (37) 

n 

For the B meson Co in the fit above is equal to ((i)\H\B)) 2 /2Eq, with a relativistic 
normalisation of the states. On the right-hand-side of Figure|22]we illustrate the operator 
in which the B meson is destroyed by the axial vector current (the W decay to leptons is 
handled analytically and separately from the lattice QCD calculation). Then 

!<0|tft( T )J^(0)|0) =« fl+(T)J Af ,(0) »=Y,D n e- E - aT . (38) 

n 

We again have a product of M^ 1 factors to average over configurations, since J ^ con- 
tains the same kind of product of ip and ip fields as H. The same energies appear in this 
fit (indeed the correlators should be fit simultaneously to ensure that this is true) but the 
amplitudes are different. Now D = ({0\H\B})({B\ J Afi |0))/2£ , . The second factor is 
the matrix element that we want and we can extract this as 2EqDq/ J (2EqCq). If the 
temporal axial current is used for a B meson at rest then /s^/msa 3 / 2 = y/2/CoDo. 

The current J a that we use on the lattice needs to be well-matched to the continuum 
current. The leading term in the lattice version of the current will be just be the obvious 
transcription of the continuum current, — ^(,7^73 V'u- However this current has dis- 
cretisation errors at 0(a s a) and these can be improved by adding higher order operators 
to cancel the errors in the same way as for the action earlier. Again we can use perturba- 
tion theory to calculate the coefficients of these correction terms. If we use NRQCD for 
the b quark in the current then there are also relativistic corrections that can be applied to 
make the current a more accurate version of the continuum current. In fact the relativistic 
correction operators and the discretisation correction operators are the same and simply 
pick up perturbative coefficients that are functions of mqa (Morningstar and Shigemitsu 
1998). 

We also have to renormalise the matrix elements from the lattice to an appropriate 
continuum renormalisation scheme such as MS. This can be done in perturbation theory 
as, again, it takes account of gluons with momenta above the lattice cut-off. Such cal- 
culations have only been done through 0(a s ) so far and this means that the final quoted 
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result has errors at 0(a s ) 2 (Morningstar and Shigemitsu 1998). Now that the system- 
atic error from working in the quenched approximation has been overcome, this is often 
the largest source of error and much more work must be done in future to reduce this 
error. Methods for renormalisation and matching that use direct numerical methods on 
the lattice (often called nonperturbative) are also being explored by many people. Matrix 
elements of conserved currents do not need renormalisation and this explains why, for 
example, the calculations of /„. and fjc described above are so accurate. 

|— P — | — i — | — i — | — i — | — i — | — i — | — i — | — i — | — i — | — i — I 
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Figure 23. The ratio of Jb s \J m b7/ '/sV m B from unquenched lattice calculations on 
the MILC configurations at two values of the lattice spacing ( Gray et al. 2005). The lines 
represent fits using chiral perturbation theory to various combinations of valence and 
sea light quark masses. The pink curve is the full QCD curve that extrapolates to the 
physical answer. 

The B meson decay constant is of interest, both because it sets the rate of B leptonic 
decay but also because it appears in the mixing rate of neutral B mesons. The mixing rate 
is parameterised by J%Bb where Bb is the 'bag constant'. J%Bb is being calculated 
in lattice QCD but it is harder than the calculation of fs alone. We believe that the Bb 
factor is fairly benign with a value around 1 . The calculation of Jb then provides a good 
indicator of the size of mixing effects (for the orange and yellow rings in Figurel2Tl until 
JbBb is known accurately. 

Figure [23] shows the current results by the HPQCD collaboration, using the MILC 
unquenched configurations and the NRQCD formalism for the valence b quarks and asq- 
tad improved staggered formalism for the valence light quarks (Gray et al. 2005, Wingate 
2004). What is plotted is the ratio of fs s yfn^ s I f B\frn~B in which the overall renormal- 
isation constant for the lattice Ja cancels giving an accurate result. What is interesting 
here is the approach to the light quark mass limit in which m u u takes its physical value. 
Chiral perturbation theory expects fairly strong logarithmic dependence but earlier results 
worked at such heavy rn u /d that they were not in the region in which chiral perturbation 
theory was valid. It is clear that the new results now have light m u /d and an accurate 
result for this ratio will be possible. 
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Figure 24. The calculation of a 3 -point function for B — > tt semileptonic decay on the 
lattice. 




Figure 25. Results for the form factors for B — > n semileptonic decay from unquenched 
lattice QCD (Shigemitsu et al. 2005). Results are shown for a variety of light quarks, 
given in terms of the mass of the tt made from these quarks. The bursts show the (small) 
extrapolation to the real tt mass. 



The determination of semileptonic decay rates requires the calculation of a 3-point 
function on the lattice. This is illustrated in Figure [24] for B — > tt. We now have two 
hadron operators at and T for different hadrons and an intermediate current operator 
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Jy^ or J a which causes the quark flavour change and the emission of a W. Now 

^(0\H*(T)J(t)H'(0)\0) =« H\T)J(t)H'(0) »= D n e' E - a{T - e - ] e~ E ' at . 

n 

(39) 

The amplitudes, D n , are now related to < 0\H\B X B\J\ix X n\H'\0 > and the 
matrix element that we want, ie < B\J\tt > can be extracted by simultaneously fitting 
the relevant 2-point functions to determine the other factors in Do- The matrix element is 
now a function of the 4-momentum transfer between B and ir, q 2 . It is expressed in terms 
of two different form factors, f+(q 2 ) and fo(q 2 ), with different momentum-dependent 
prefactors. f+(q 2 ) is the form factor that translates directly into the rate of semileptonic 
decay, since /q ends up in this rate multiplied by the mass of the lepton into which the 
W decays, which is very small. 
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Figure 26. Results for the form factors for D semileptonic decay to tt and K. The 
experimental points shown at q 2 — use the experimental decay rate at that point and 
the current PDG results for the appropriate CKM elements (Aubin et al. 2005). 

Once again we must match the lattice current to a continuum current and renormalise 
to obtain results in a continuum renormalisation scheme. This limits the accuracy with 
which these calculations can now be done and more work is required to improve this 
situation. Figure |25] and Figure [26] show current results obtained by the HPQCD and 
FNAL/MILC collaborations for the formfactors for B — > tt decay and D — > tt/K decay 
respectively. The B results use NRQCD b quarks (Shigemitsu et al. 2005); the D results 
use the Fermilab formalism for the c quark (Aubin et al. 2005). The usefulness of the D 
results is to compare to imminent experimental results from the CLEO-c collaboration 
that will provide a check of lattice methods and systematic errors for confidence in our 
precision B results (Shipsey 2005). 



Lattice QCD 

4 Conclusions 



39 



Lattice QCD has come a long way from the original calculations of the 1970s. The origi- 
nal idea that we could solve a simple discretisation of QCD numerically by 'brute force' 
has been replaced by a more sophisticated approach. Improved discretisation of both the 
gluon action and the quark action has led to the possibility of performing realistic simu- 
lations of QCD on current computers. Indeed this has been done, and I hope that I have 
conveyed something of the excitement of seeing accurate lattice calculations reproduce 
well-known experimental numbers for the first time. These first calculations are of the 
masses of 'gold-plated' hadrons, those for which lattice QCD must be able to get the 
right answer if it is to be trusted at all. Leading on from this we have been able to make 
the first accurate lattice QCD prediction of the mass of a new meson, the B c . 

It is now important to beat down the sources of systematic error in the lattice calcula- 
tion of decay matrix elements for B and D physics to obtain results that can be combined 
with experiment to give an accurate determination of elements of the CKM matrix. The 
timescale for this programme is the next two years and, as well as lattice QCD calcu- 
lations, it requires 0{a 2 s ) calculations in perturbation theory to renormalise the lattice 
results to numbers appropriate to the continuous real world. On a longer timescale (say 
five years) studies of more complicated baryonic matrix elements will be undertaken and 
progress will be made in understanding to what extent accurate lattice calculations can 
be done of some of the more interesting, but not gold-plated, particles in the hadron 
spectrum. 
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